Adaptive potential Merge genetic data
Bioshifts group
1 Load libraries and source code
# devtools::install_github('EcologicalTraitData/traitdataform')
rm(list=ls())
list.of.packages <- c("ggplot2","data.table","dplyr","tidyr","parallel","bdc","taxadb","traitdataform","pbapply","tidyverse","readxl","lme4","coefplot","sjPlot","sjmisc","effects","rgdal","maptools","rgeos","terra","MuMIn","rnaturalearthdata","lsmeans","GGally","tidyterra","httr","purrr","rlist","usethis")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)
sapply(list.of.packages, require, character.only = TRUE)## Warning: package 'ggplot2' was built under R version 4.3.2## Warning: package 'terra' was built under R version 4.3.2## Warning: package 'tidyterra' was built under R version 4.3.2## ggplot2 data.table dplyr tidyr
## TRUE TRUE TRUE TRUE
## parallel bdc taxadb traitdataform
## TRUE TRUE TRUE TRUE
## pbapply tidyverse readxl lme4
## TRUE TRUE TRUE TRUE
## coefplot sjPlot sjmisc effects
## TRUE TRUE TRUE TRUE
## rgdal maptools rgeos terra
## TRUE TRUE TRUE TRUE
## MuMIn rnaturalearthdata lsmeans GGally
## TRUE TRUE TRUE TRUE
## tidyterra httr purrr rlist
## TRUE TRUE TRUE TRUE
## usethis
## TRUE2 source functions
# Code to find accepted species names
# Before you can download the source code from github, make sure you have a personal github token
# run this:
# usethis::edit_r_environ()
# if you have no toke, create one following:
# 1. Generate on GitHub your personal token
# 1.1 Go to GitHub
# 2.1 In the right corner go to "Settings"
# 2.2 Then in the left part go to "Developer setting"
# 2.3 Select the option "Personal access tokens"
# 2.4 Select the option "Generate new token"
# 2.5 Copy your personal token
# run this to add the token to your .Renviron
# usethis::edit_r_environ()
# write GITHUB_PAT=YOUR_TOKEN
# Sys.getenv("GITHUB_PAT")
source(here::here("R/fetchGHdata.R"))
fetchGHdata(gh_account = "Bioshifts",
repo = "bioshifts_v1_v2",
path = "R/Source_code/Clean_names.R",
output = here::here("R/Clean_names.R"))## Loading required package: jsonlite##
## Attaching package: 'jsonlite'## The following object is masked from 'package:purrr':
##
## flatten## [1] 0fetchGHdata(gh_account = "Bioshifts",
repo = "bioshifts_v1_v2",
path = "R/Source_code/Find_Sci_Names.R",
output = here::here("R/Find_Sci_Names.R"))## [1] 0source(here::here("R/Clean_names.R"))
source(here::here("R/Find_Sci_Names.R"))
# Malin's transformation for moving values slightly inward.
transform01 <- function(x) (x * (length(x) - 1) + 0.5) / (length(x))
# De Kort transformation
deKort_trans <- function(p){
p <- scale(p)
p <- (p - min(p, na.rm = T))/(max(p, na.rm = T)-min(p, na.rm = T))
return(p)
}
# Malin's transformation before a logit transformation
logit_trans <- function(p){
p <- transform01(p)
p <- log(p/(1-p))
return(p)
}3 Load species list
splist <- read.csv(here::here("Data/splist.csv"), header = T)
# remove duplicated sp_names
splist <- splist %>%
dplyr::filter(!duplicated(scientificName))%>%
mutate(Kingdom=kingdom,
Phylum=phylum,
Class=class,
Order=order,
Family=family)%>%
dplyr::select(scientificName,Kingdom,Phylum,Class,Order,Family,db)
splist$Genus <- sapply(splist$scientificName, function(x){
strsplit(x," ")[[1]][1]
})4 Load bioshifts
4.1 Geo data
if(!file.exists(here::here("data/bioshiftsv3_metadata.csv"))){
# # From Google drive
# drive_path <- readWindowsShortcut("G:/My Drive.lnk")$pathname
# bioshifts_path <- readWindowsShortcut(paste0(drive_path,"\\BIOSHIFTS Working Group.lnk"))$pathname
# shp_path <- here::here(bioshifts_path,("DATA/BIOSHIFTS_v3/CleanDatabasev3/ShapefilesBioShiftsv3"))
# all_shps <- list.files(shp_path,pattern = "dbf", full.names = TRUE)
# From local folder
all_shps <- list.files("C:/Users/brunn/ShadowDrive/CreateGeodatabaseBioShifts/Data/ShapefilesBioShiftsv3",
pattern = "shp", full.names = TRUE)
biov1_meta <- pblapply(all_shps, terra::vect)
biov1_meta <- pblapply(biov1_meta, data.frame)
biov1_meta <- data.table::rbindlist(biov1_meta)
write.csv(biov1_meta, here::here("data/bioshiftsv3_metadata.csv"),row.names = FALSE)
} else {
biov1_meta <- read.csv(here::here("data/bioshiftsv3_metadata.csv"))
}4.2 Load v1
biov1 <- read.csv(here::here("Data/biov1_fixednames.csv"), header = T)
# Fix references in biov1
biov1$sp_name_std_v1 <- gsub("_"," ",biov1$sp_name_std_v1)
# Select columns of interest
biov1 <- biov1 %>%
filter(!Type %in% c("BAT","LON")) %>%
dplyr::select(ID,Article_ID,Study_ID,
Type,Param,Trend,SHIFT,UNIT,DUR,
v.lat.mean,v.ele.mean,
START,END,Sampling,Uncertainty_Distribution,Uncertainty_Parameter,Nperiodes,
N,Grain_size,Data,ID.area,
Phylum,Class,Order,Family,Genus,sp_name_std_v1,
group,ECO,Hemisphere) %>% # select columns
mutate(
Type = case_when(
Type=="HOR" ~ "LAT",
TRUE ~ as.character(Type)),
Data = case_when(
Data=="occurence-based" ~ "occurrence-based",
TRUE ~ as.character(Data)),
spp = sp_name_std_v1,
SHIFT_abs = abs(SHIFT),
velocity = ifelse(Type == "LAT", v.lat.mean, v.ele.mean),
vel_sign = ifelse(velocity>0,"pos","neg"))
# add geospatial data
biov1_meta$ID <- biov1_meta$Name
biov1 <- biov1 %>%
filter(ID %in% biov1_meta$ID)
# merge
biov1 <- merge(biov1,biov1_meta,by="ID",all.y = TRUE) # all.y to keep only studies with geospatial data
# Get Extent (=Lat/Ele extent)
biov1$Extent <- biov1$LatExtentk
biov1$Extent[which(biov1$Type=="ELE")] <- biov1$EleExtentm[which(biov1$Type=="ELE")]
# Number of temporal units
biov1$NtempUnits <- biov1$Nperiodes
these <- which(biov1$NtempUnits > biov1$DUR)
biov1$NtempUnits[these] <- biov1$DUR[these]
# log NtempUnits
biov1$LogNtempUnits <- log(biov1$NtempUnits)
## group irregular sampling with multiple sampling
biov1$Sampling = ifelse(biov1$Sampling %in% c("IRR","MULTIPLE(continuous)"),"MULT", biov1$Sampling)
## Grain
biov1$Grain <- factor(biov1$Grain_size, levels = c("small", "moderate", "large","very_large"))
biov1$ContinuousGrain <-
as.numeric(
as.character(
ifelse(
biov1$Grain == "small", 1,
ifelse(biov1$Grain == "moderate",
2, #10,
ifelse(biov1$Grain == "large",
3, #100,
4))))) #we have a few very large in the full database
#harmonize column names
biov1$Quality <- ifelse(biov1$Uncertainty_Distribution %in% c("RESAMPLING"),"BALANCED",
ifelse(biov1$Uncertainty_Distribution %in% c("MODEL","MODEL+RESAMPLING(same)","RESAMPLING(same)+MODEL","RESAMPLING+MODEL"),"BALANCED",
ifelse(biov1$Uncertainty_Distribution %in% c("DETECTABILITY","RESAMPLING(same)+DETECTABILITY"),"BALANCED",
ifelse(biov1$Uncertainty_Distribution %in% c("RESAMPLING(same)"),"RESURVEYED",biov1$Uncertainty_Distribution))))
table(biov1$Quality) ##
## BALANCED RAW RESURVEYED
## 14050 9497 6522biov1$PrAb <- ifelse(biov1$Data %in% c("occurence-based", "occurrence-based"),"occurrence","abundance")
table(biov1$PrAb) ##
## abundance occurrence
## 9544 20571biov1 <- biov1 %>%
dplyr::filter(!is.na(sp_name_std_v1))
all(biov1$sp_name_std_v1 %in% splist$scientificName)## [1] TRUE# from continuous to categorical
q1=quantile(biov1$START,probs=c(0,0.25,0.5,0.75,1))
biov1$StartF=cut(biov1$START,breaks=q1,include.lowest=T)
q1=quantile(biov1$Extent ,probs=c(0,0.25,0.5,0.75,1))
biov1$ExtentF=cut(biov1$Extent,breaks=q1,include.lowest=T)
q1=quantile(biov1$ID.area,probs=c(0,0.25,0.5,0.75,1))
biov1$AreaF=cut(biov1$ID.area,breaks=q1,include.lowest=T)
q1=quantile(biov1$N,probs=c(0,0.25,0.5,0.75,1))
biov1$NtaxaF=cut(biov1$N,breaks=q1,include.lowest=T)
# add ID to obs
biov1$obs_ID <- paste0("S",1:nrow(biov1))
summary(biov1$velocity)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## -5.8328 0.5415 2.0271 2.3071 2.8319 13.1259 1summary(biov1$velocity[which(biov1$vel_sign=="pos")])## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.009396 0.541507 2.055833 2.434770 2.831932 13.1258854.3 Classify species
# Class species as Terrestrial, Marine or Freshwater
Terv1 <- unique(biov1$sp_name_std_v1[which(biov1$ECO == "T")])
Terrestrials = unique(c(Terv1))
Mar1 <- unique(biov1$sp_name_std_v1[which(biov1$ECO == "M")])
Marine = unique(c(Mar1))
# Freshwater fish
FFishv1 <- unique(biov1$sp_name_std_v1[(biov1$Class == "Actinopterygii" | biov1$Class == "Cephalaspidomorphi") & biov1$ECO == "T"])
FFish = unique(c(FFishv1))
# Marine fish
MFishv1 <- unique(biov1$sp_name_std_v1[biov1$Class == "Actinopterygii" & biov1$ECO == "M"])
MFish = unique(c(MFishv1))
splist$ECO = NA
splist$ECO[which(splist$scientificName %in% Terrestrials)] <- "T"
splist$ECO[which(splist$scientificName %in% Marine)] <- "M"
splist$ECO[which(splist$scientificName %in% MFish)] <- "M"
biov1$ECO = NA
biov1$ECO[which(biov1$sp_name_std_v1 %in% Terrestrials)] <- "T"
biov1$ECO[which(biov1$sp_name_std_v1 %in% Marine)] <- "M"
biov1$ECO[which(biov1$sp_name_std_v1 %in% MFish)] <- "M"
splist$Group = NA
splist$Group[which(splist$Class == "Phaeophyceae")] <- "Chromista"
splist$Kingdom[which(splist$Class == "Phaeophyceae")] <- "Chromista"
splist$Group[which(splist$Phylum == "Rhodophyta")] <- "Seaweed"
splist$Kingdom[which(splist$Phylum == "Rhodophyta")] <- "Plantae"
splist$Group[which(splist$Family == "Elminiidae")] <- "Barnacles"
splist$Group[which(splist$Kingdom == "Bacteria")] <- "Bacteria"
splist$Group[which(splist$Class == "Holothuroidea")] <- "Sea cucumber"
splist$Group[which(splist$Class == "Aves")] <- "Bird"
splist$Group[which(splist$Class == "Insecta")] <- "Insect"
splist$Group[which(splist$Class == "Mammalia")] <- "Mammal"
splist$Group[which(splist$Class == "Arachnida")] <- "Spider"
splist$Kingdom[which(splist$Kingdom == "Viridiplantae")] <- "Plantae"
splist$Kingdom[which(splist$Phylum == "Tracheophyta")] <- "Plantae"
splist$Group[which(splist$Kingdom == "Plantae")] <- "Plant"
splist$Group[which(splist$Class == "Hydrozoa")] <- "Hydrozoa"
splist$Group[which(splist$Class == "Anthozoa")] <- "Sea anemones and corals"
splist$Group[which(splist$Class == "Polychaeta")] <- "Polychaetes"
splist$Group[which(splist$Phylum == "Mollusca")] <- "Molluscs"
splist$Group[which(splist$Class == "Malacostraca")] <- "Crustacean"
splist$Group[which(splist$Class == "Hexanauplia")] <- "Crustacean"
splist$Group[which(splist$Class == "Maxillopoda")] <- "Crustacean"
splist$Group[which(splist$Class == "Ostracoda")] <- "Crustacean"
splist$Group[which(splist$Class == "Branchiopoda")] <- "Crustacean"
splist$Group[which(splist$Class == "Asteroidea")] <- "Starfish"
splist$Group[which(splist$Class == "Ascidiacea")] <- "Ascidians tunicates and sea squirts"
splist$Class[which(splist$Class == "Actinopteri")] <- "Actinopterygii"
splist$Group[which(splist$Class == "Actinopterygii")] <- "Fish"
splist$Group[which(splist$Class == "Elasmobranchii")] <- "Fish"
splist$Group[which(splist$Order == "Perciformes")] <- "Fish"
splist$Group[which(splist$Class == "Chondrichthyes")] <- "Fish"
splist$Group[which(splist$Class == "Holocephali")] <- "Fish"
splist$Group[which(splist$Class == "Cephalaspidomorphi")] <- "Fish"
splist$Group[which(splist$Class == "Echinoidea")] <- "Sea urchin"
splist$Group[which(splist$Class == "Crinoidea")] <- "Crinoid"
splist$Group[which(splist$Class == "Holothuroidea")] <- "Sea cucumber"
splist$Group[which(splist$Class == "Reptilia")] <- "Reptile"
splist$Group[which(splist$Order == "Squamata")] <- "Reptile"
splist$Group[which(splist$Class == "Ophiuroidea")] <- "Brittle stars"
splist$Group[which(splist$Class == "Chilopoda")] <- "Centipedes"
splist$Group[which(splist$Class == "Diplopoda")] <- "Millipedes"
splist$Group[which(splist$Class == "Amphibia")] <- "Amphibian"
splist$Group[which(splist$Kingdom == "Fungi")] <- "Fungi"
splist$Group[which(splist$Order == "Balanomorpha")] <- "Barnacles"
splist$Group[which(splist$Phylum == "Nematoda")] <- "Nematodes"
splist$Group[which(splist$Class == "Myxini")] <- "Hagfish"
splist$Group[which(splist$Kingdom == "Chromista")] <- "Chromista"
splist$Family[which(splist$Genus == "Dendrocopus")] <- "Picidae"
######################################
biov1 <- merge(biov1[,-which(names(biov1) %in% c("Phylum","Class","Order","Family","Genus","Group","ECO"))],
splist[,c("Kingdom","Phylum","Class","Order","Family","Genus","Group","ECO","scientificName")],
by.x = "sp_name_std_v1", by.y = "scientificName",
all.x = T)4.4 Filter LAT / ELE shifts
Use only latitudinal and elevational shifts
# v1
biov1 <- biov1 %>%
dplyr::filter(Type %in% c("ELE","LAT")) # Use LAT ELE shifts
# splist
sps <- unique(biov1$sp_name)
splist <- splist %>% dplyr::filter(scientificName %in% sps)4.5 Remove species identified to the genus level or cf.
if(any(grep("sp[.]",biov1$sp_reported_name_v1))){
biov1 <- biov1 %>% dplyr::filter(!grepl("sp[.]",sp_reported_name_v1))
}
if(any(grep("sp[.]",biov1$sp_name_std_v1))){
biov1 <- biov1 %>% dplyr::filter(!grepl("sp[.]",sp_name_std_v1))
}
if(any(grep("cf[.]",biov1$sp_reported_name_v1))){
biov1 <- biov1 %>% dplyr::filter(!grepl("cf[.]",sp_reported_name_v1))
}
if(any(grep("cf[.]",biov1$sp_name_std_v1))){
biov1 <- biov1 %>% dplyr::filter(!grepl("cf[.]",sp_name_std_v1))
}
#remove freshwater fishes
biov1 <- biov1[-which((biov1$Class == "Actinopterygii" | biov1$Group == "Fish") & biov1$ECO=="T"),]
#remove marine birds
biov1 <- biov1[-which(biov1$Class == "Aves" & biov1$ECO=="M"),]
if(any(grep("sp[.]",splist$scientificName))){
splist <- splist %>% dplyr::filter(!grepl("sp[.]",scientificName))
}
if(any(grep("sp[.]",splist$scientificName))){
splist <- splist %>% dplyr::filter(!grepl("sp[.]",scientificName))
}
if(any(grep("cf[.]",splist$scientificName))){
splist <- splist %>% dplyr::filter(!grepl("cf[.]",scientificName))
}
if(any(grep("cf[.]",splist$scientificName))){
splist <- splist %>% dplyr::filter(!grepl("cf[.]",scientificName))
}4.6 Reclassify methodological variables
unique(biov1$Data)## [1] "abundance-based" "occurrence-based"biov1$Data[biov1$Data=="occurence-based"] <- "occurrence-based"
biov1$Data <- factor(biov1$Data, levels = unique(biov1$Data))
table(biov1$Data)##
## abundance-based occurrence-based
## 9421 20393unique(biov1$Sampling)## [1] "TWO" "CONTINUOUS" "MULT"biov1$Sampling = ifelse(biov1$Sampling %in% c("IRR","MULTIPLE(continuous)"),"MULTIPLE", biov1$Sampling)
biov1$Sampling <- ordered(biov1$Sampling,
levels = c("TWO","MULTIPLE","CONTINUOUS"))
table(biov1$Sampling)##
## TWO MULTIPLE CONTINUOUS
## 25466 0 3105unique(biov1$Grain_size)## [1] "small" "large" "moderate" "very_large" NAbiov1$Grain_size <- ifelse(biov1$Grain_size %in% c("large","very_large"),"large",biov1$Grain_size)
biov1$Grain_size <- ordered(biov1$Grain_size,
levels = c("small","moderate","large"))
table(biov1$Grain_size)##
## small moderate large
## 9740 10657 9416unique(biov1$Uncertainty_Distribution)## [1] "RESAMPLING(same)" "RAW"
## [3] "RESAMPLING" "MODEL"
## [5] "RESAMPLING+MODEL" "MODEL+RESAMPLING(same)"
## [7] "RESAMPLING(same)+DETECTABILITY" "RESAMPLING(same)+MODEL"
## [9] "DETECTABILITY"biov1$Uncertainty_Distribution <- ifelse(biov1$Uncertainty_Distribution %in% c("RESAMPLING","RESAMPLING(same)"),"RESAMPLING",
ifelse(biov1$Uncertainty_Distribution %in% c("MODEL","MODEL+RESAMPLING(same)","RESAMPLING+MODEL"),"MODEL",
ifelse(biov1$Uncertainty_Distribution %in% c("DETECTABILITY","RESAMPLING(same)+DETECTABILITY"),"DETECTABILITY",
biov1$Uncertainty_Distribution)))
biov1$Uncertainty_Distribution <- ifelse(biov1$Uncertainty_Distribution == "RAW","OPPORTUNISTIC","PROCESSED")
table(biov1$Uncertainty_Distribution)##
## OPPORTUNISTIC PROCESSED
## 9449 20365#transform study area
biov1$ID.area <- log(biov1$ID.area)4.7 Get centroids of study areas
This will be used to merge genetic data with bioshifts based on distance of observations
# # The input file geodatabase
# fgdb <- "C:/Users/brunn/nextCloud/bioshifts_v1_v2/v1/Study_Areas_v1/Study_Areas.gdb"
#
# # List all feature classes in a file geodatabase
# fc_list <- ogrListLayers(fgdb)
#
# # Get centroids
# cl <- makeCluster(detectCores()-1)
# clusterExport(cl, c("readOGR", "gCentroid", "fgdb"))
#
# centroids <- pblapply(fc_list, function(x){
# fc <- readOGR(dsn=fgdb,layer=x)
# c <- data.frame(gCentroid(fc))
# names(c) <- c("centroid.x", "centroid.y")
#
# geomet <- data.frame(terra::geom(terra::vect(fc)))
#
# max_lat <- order(geomet[,"y"], decreasing = T)[1]
# max_lat <- geomet[max_lat,c("x","y")]
# names(max_lat) <- c("max_lat.x", "max_lat.y")
#
# min_lat <- order(geomet[,"y"])[1]
# min_lat <- data.frame(geomet[min_lat,c("x","y")])
# names(min_lat) <- c("min_lat.x", "min_lat.y")
#
# cbind(c, max_lat, min_lat)
# }, cl = cl)
#
# stopCluster(cl)
#
# centroids <- do.call("rbind",centroids)
# centroids$ID <- fc_list
#
# write.csv(centroids, "Data/centroids_study_areas.csv", row.names = FALSE)
centroids <- read.csv(here::here("Data/centroids_study_areas.csv"))
biov1 <- merge(biov1, centroids, by = "ID")5 Corrected shifts
Just load it in
fetchGHdata(gh_account = "Bioshifts",
repo = "MethodologicalAdjustment",
path = "outputs/biov1_method_corrected_shifts_study_level.csv",
output = here::here("Data/biov1_method_corrected_shifts_study_level.csv"))## [1] 0biov1_corr <- read.csv(here::here("Data/biov1_method_corrected_shifts_study_level.csv"))
# add corrected shifts
biov1 <- merge(biov1,
biov1_corr %>%
select(c("Article_ID","Study_ID","Class","Type","SLDiff1")),
by = c("Article_ID","Study_ID","Class","Type"),
all.x = TRUE)
biov1$SHIFT_cor <- abs(biov1$SHIFT) - biov1$SLDiff1
# change sign for negative shifts
biov1$neg_shifts <- ifelse(sign(biov1$SHIFT) == -1 & biov1$SHIFT != 0, 1, 0)
biov1$SHIFT_cor[which(biov1$neg_shifts == 1)] <- biov1$SHIFT_cor[which(biov1$neg_shifts == 1)] * -16 Direction of shift
Identify direction of shift
biov1 <- biov1 %>%
mutate(shift_sign = ifelse(SHIFT>0,"pos","neg"),
shift_vel_sign = paste0(shift_sign,vel_sign),
shiftC_sign = ifelse(SHIFT_cor>0,"pos","neg"),
shiftC_vel_sign = paste0(shiftC_sign,vel_sign)
) %>%
filter(!is.na(velocity))
ggplot(biov1, aes(x=shift_vel_sign))+
geom_bar()+
coord_flip()+
facet_wrap(Type~Param)
ggplot(biov1, aes(x=shiftC_vel_sign))+
geom_bar()+
coord_flip()+
facet_wrap(Type~Param)
summary(biov1$velocity)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -5.8328 0.5415 2.0271 2.2860 2.7543 13.1259summary(biov1$velocity[which(biov1$shift_vel_sign=="pospos")])## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.009396 0.541507 2.073300 2.470806 2.831932 13.125880summary(biov1$velocity[which(biov1$shift_vel_sign=="negneg")])## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -3.79347 -2.30487 -1.78004 -1.59872 -0.58152 -0.029457 Calculate lags raw
# calculate lags
# positive values are a lag (range shift lower smaller than expected or in opposite directions) and negative values are range shift larger than expected
biov1$lag <- biov1$velocity-biov1$SHIFT
position = which(biov1$vel_sign == "neg")
biov1$lag[position] <- biov1$lag[position] * -1 # Any negative velocity means the sign of lag has to shift.
# Lag 2
# When shift is in opposite sign of velocity, lag = abs(velocity)
biov1$lag2 = biov1$lag
position <- which(biov1$shift_vel_sign == "posneg" | biov1$shift_vel_sign == "negpos")
biov1$lag2[position] <- abs(biov1$velocity[position])
# Lag 3
# When shift is > velocity, lag = 0
biov1$lag3 = biov1$lag2
position <- which((biov1$shift_sign == "pos" & biov1$SHIFT > biov1$velocity) |
(biov1$shift_sign == "neg" & biov1$SHIFT < biov1$velocity))
biov1$lag3[position] <- 0
{
par(mfrow=c(2,2))
plot(lag~velocity, biov1)
plot(lag2~velocity, biov1)
plot(lag3~velocity, biov1)
}
ggplot(biov1, aes(x = Param, y=lag, fill = Param, color = Param))+
geom_boxplot(alpha = .5, outlier.shape = NA)+
scale_y_continuous(limits = quantile(biov1$lag2, c(0.1, 0.9), na.rm = T))+
facet_wrap(.~Type, scales = "free")## Warning: Removed 8756 rows containing non-finite values (`stat_boxplot()`).ggplot(biov1, aes(x = Param, y=lag2, fill = Param, color = Param))+
geom_boxplot(alpha = .5, outlier.shape = NA)+
scale_y_continuous(limits = quantile(biov1$lag2, c(0.1, 0.9), na.rm = T))+
facet_wrap(.~Type, scales = "free")## Warning: Removed 5850 rows containing non-finite values (`stat_boxplot()`).ggplot(biov1, aes(x = Param, y=log1p(lag3), fill = Param, color = Param))+
geom_boxplot(alpha = .5, outlier.shape = NA)+
# scale_y_continuous(limits = quantile(log1p(biov1$lag3), c(0.1, 0.9), na.rm = T))+
facet_wrap(.~Type, scales = "free")
8 Calculate lags corrected
# calculate lagCs
# positive values are a lagC (range shift lower smaller than expected or in opposite directions) and negative values are range shift larger than expected
biov1$lagC <- biov1$velocity-biov1$SHIFT_cor
position = which(biov1$vel_sign == "neg")
biov1$lagC[position] <- biov1$lagC[position] * -1 # Any negative velocity means the sign of lagC has to shift.
# Lag 2
# When shift is in opposite sign of velocity, lagC == velocity
biov1$lagC2 = biov1$lagC
position <- which(biov1$shiftC_vel_sign == "posneg" | biov1$shiftC_vel_sign == "negpos")
biov1$lagC2[position] <- abs(biov1$velocity[position])
# Lag 3
# When shift is > velocity, lagC = 0
biov1$lagC3 = biov1$lagC2
position <- which((biov1$shiftC_sign == "pos" & biov1$SHIFT_cor > biov1$velocity) |
(biov1$shiftC_sign == "neg" & biov1$SHIFT_cor < biov1$velocity))
biov1$lagC3[position] <- 0
{
par(mfrow=c(2,2))
plot(lagC~velocity, biov1)
plot(lagC2~velocity, biov1)
plot(lagC3~velocity, biov1)
}
ggplot(biov1, aes(x = Param, y=lagC, fill = Param, color = Param))+
geom_boxplot(alpha = .5, outlier.shape = NA)+
scale_y_continuous(limits = quantile(biov1$lagC2, c(0.1, 0.9), na.rm = T))+
facet_wrap(.~Type, scales = "free")## Warning: Removed 10339 rows containing non-finite values (`stat_boxplot()`).ggplot(biov1, aes(x = Param, y=lagC2, fill = Param, color = Param))+
geom_boxplot(alpha = .5, outlier.shape = NA)+
scale_y_continuous(limits = quantile(biov1$lagC2, c(0.1, 0.9), na.rm = T))+
facet_wrap(.~Type, scales = "free")## Warning: Removed 6147 rows containing non-finite values (`stat_boxplot()`).ggplot(biov1, aes(x = Param, y=log1p(lagC3), fill = Param, color = Param))+
geom_boxplot(alpha = .5, outlier.shape = NA)+
# scale_y_continuous(limits = quantile(log1p(biov1$lagC3), c(0.1, 0.9), na.rm = T))+
facet_wrap(.~Type, scales = "free")## Warning: Removed 229 rows containing non-finite values (`stat_boxplot()`).
9 Corrected vs raw shifts
plot(SHIFT_cor~SHIFT,biov1)
abline(a=0,b=1,col=2)
lm0=lm(SHIFT_cor~SHIFT,biov1)
summary(lm0)##
## Call:
## lm(formula = SHIFT_cor ~ SHIFT, data = biov1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.4752 -1.0085 0.1219 0.8303 3.4359
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.112782 0.007375 15.29 <2e-16 ***
## SHIFT 1.020616 0.001328 768.28 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.241 on 29582 degrees of freedom
## (229 observations deleted due to missingness)
## Multiple R-squared: 0.9523, Adjusted R-squared: 0.9523
## F-statistic: 5.903e+05 on 1 and 29582 DF, p-value: < 2.2e-16hist(biov1$SHIFT,col=rgb(1,0,0,0.5))
hist(biov1$SHIFT_cor,col=rgb(0,0,1,0.5),add=T)
#so by using residuals we decrease the variation in the raw range shift obs, it's like we work on a more homogenous variables as all the variations due to methods variation has been substracted to the shift.
x1=tapply(biov1$SHIFT,biov1$Group,mean)
x2=tapply(biov1$SHIFT_cor,biov1$Group,mean)#lower mean value than true obs
plot(x1~x2, xlab="Mean shift per Group (raw)", ylab="Mean shift per Group (corrected)")
lm0=lm(x1~x2,biov1)
summary(lm0) #it changes many things##
## Call:
## lm(formula = x1 ~ x2, data = biov1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.84859 -0.10546 0.02355 0.35682 0.47638
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.398251 0.118757 -3.353 0.00574 **
## x2 1.014730 0.009711 104.497 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4436 on 12 degrees of freedom
## (10 observations deleted due to missingness)
## Multiple R-squared: 0.9989, Adjusted R-squared: 0.9988
## F-statistic: 1.092e+04 on 1 and 12 DF, p-value: < 2.2e-16x1=tapply(biov1$SHIFT,biov1$Group,var)
x2=tapply(biov1$SHIFT_cor,biov1$Group,var)#lower variance value than in true obs
plot(x1~x2, xlab="Variance shift per Group (raw)", ylab="Variance shift per Group (corrected)")
lm0=lm(x1~x2,biov1)
summary(lm0) #but high positive correlation meaning that relative variation is conserved##
## Call:
## lm(formula = x1 ~ x2, data = biov1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.5104 -0.9744 0.5092 1.8051 5.9561
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.85821 1.70831 -1.673 0.125
## x2 1.02278 0.03198 31.980 2.1e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.734 on 10 degrees of freedom
## (12 observations deleted due to missingness)
## Multiple R-squared: 0.9903, Adjusted R-squared: 0.9893
## F-statistic: 1023 on 1 and 10 DF, p-value: 2.104e-11x1=tapply(biov1$SHIFT,biov1$sp_name_std_v1,mean)
x2=tapply(biov1$SHIFT_cor,biov1$sp_name_std_v1,mean)#lower mean value than in true obs
plot(x1~x2, xlab="Mean shift per species (raw)", ylab="Mean shift per species (corrected)")
lm0=lm(x1~x2,biov1)
summary(lm0) #but at the species level, we observe a high positive relationship so the corrected shifts did not change the pattern of range shift=> good##
## Call:
## lm(formula = x1 ~ x2, data = biov1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.0964 -0.6970 0.0045 0.7380 8.4140
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.055072 0.008903 -6.186 6.39e-10 ***
## x2 0.949746 0.001563 607.769 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9392 on 11938 degrees of freedom
## (127 observations deleted due to missingness)
## Multiple R-squared: 0.9687, Adjusted R-squared: 0.9687
## F-statistic: 3.694e+05 on 1 and 11938 DF, p-value: < 2.2e-16#looking at the relationship with climate velocity
#WARNING:here you need to adapt the code: if you look at the elevational range shift so you have to use EleVeloT in order to use the corresponding climate velocity
par(mfrow=c(1,2))
plot(SHIFT~velocity, xlab="SHIFT (raw)", ylab="Velocity",biov1)
lm0=lm(SHIFT~velocity+I(velocity^2),biov1)
summary(lm0)##
## Call:
## lm(formula = SHIFT ~ velocity + I(velocity^2), data = biov1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -120.782 -1.655 -0.782 1.253 145.242
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.855456 0.050371 16.983 < 2e-16 ***
## velocity 0.190516 0.029308 6.501 8.13e-11 ***
## I(velocity^2) -0.010528 0.002911 -3.616 3e-04 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.529 on 29810 degrees of freedom
## Multiple R-squared: 0.002352, Adjusted R-squared: 0.002285
## F-statistic: 35.14 on 2 and 29810 DF, p-value: 5.722e-16x1=seq(min(biov1$velocity,na.rm = T),max(biov1$velocity,na.rm = T),le=100)
p1=predict(lm0,newdata=data.frame(velocity=x1),type="response")
points(p1~x1,type="l",col=2,lwd=2)
plot(SHIFT_cor~velocity, xlab="SHIFT (corrected)", ylab="Velocity",biov1)#the pattern is changing a lot
lm0=lm(SHIFT_cor~velocity+I(velocity^2),biov1)
summary(lm0)##
## Call:
## lm(formula = SHIFT_cor ~ velocity + I(velocity^2), data = biov1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -119.895 -2.436 -0.109 1.758 144.074
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.936871 0.052087 17.987 < 2e-16 ***
## velocity 0.214470 0.030694 6.987 2.86e-12 ***
## I(velocity^2) -0.011038 0.003114 -3.545 0.000394 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.67 on 29581 degrees of freedom
## (229 observations deleted due to missingness)
## Multiple R-squared: 0.003057, Adjusted R-squared: 0.002989
## F-statistic: 45.35 on 2 and 29581 DF, p-value: < 2.2e-16x1=seq(min(biov1$velocity,na.rm = T),max(biov1$velocity,na.rm = T),le=100)
p1=predict(lm0,newdata=data.frame(velocity=x1),type="response")
points(p1~x1,type="l",col=2,lwd=2)
#we observed a relationship for the raw range shift (an unimodal relationship with higher shifts when absolute climate velocity increase).
#for the range shift corrected by method, no relationship with climatic velocity is observed
#looking at the relationship with lag
par(mfrow=c(1,2))
plot(lagC~velocity,biov1,main="corrected lag estimate")
# Lag C2
# When shift is in opposite sign of velocity, lag == velocity
plot(lagC2~velocity,biov1,main="corrected lag estimate with correction for special case")
#change more things than the above lag metrics, still the highly change are observed at extreme negative values
plot(lag2~lagC2,biov1)
lm0=lm(lag~lagC,biov1)
summary(lm0) #we observe a high positive relationship so the corrected lag did not change the pattern of the observed lag=> good##
## Call:
## lm(formula = lag ~ lagC, data = biov1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.3871 -0.9069 -0.0191 0.7850 7.9208
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.179997 0.007104 25.34 <2e-16 ***
## lagC 0.943464 0.001152 818.81 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.2 on 29582 degrees of freedom
## (229 observations deleted due to missingness)
## Multiple R-squared: 0.9577, Adjusted R-squared: 0.9577
## F-statistic: 6.705e+05 on 1 and 29582 DF, p-value: < 2.2e-16par(mfrow=c(1,2))
plot(lag2~velocity,biov1)
lm0=lm(lag2~velocity+I(velocity^2),biov1)
summary(lm0) #significant; R2=46%##
## Call:
## lm(formula = lag2 ~ velocity + I(velocity^2), data = biov1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -144.654 -0.370 1.128 2.127 6.562
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.902258 0.038778 -23.27 <2e-16 ***
## velocity 0.283622 0.022563 12.57 <2e-16 ***
## I(velocity^2) 0.053700 0.002241 23.96 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.256 on 29810 degrees of freedom
## Multiple R-squared: 0.1901, Adjusted R-squared: 0.19
## F-statistic: 3498 on 2 and 29810 DF, p-value: < 2.2e-16x1=seq(min(biov1$velocity, na.rm = T),max(biov1$velocity, na.rm = T),le=100)
p1=predict(lm0,newdata=data.frame(velocity=x1),type="response")
points(p1~x1,type="l",col=2,lwd=2) #bimodal relationship
plot(lagC2~velocity,biov1)
lm0=lm(lagC2~velocity+I(velocity^2),biov1)
summary(lm0)#significant; R2=86##
## Call:
## lm(formula = lagC2 ~ velocity + I(velocity^2), data = biov1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -143.281 -0.624 1.177 2.150 6.533
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.160582 0.038835 -29.89 <2e-16 ***
## velocity 0.235309 0.022885 10.28 <2e-16 ***
## I(velocity^2) 0.053876 0.002322 23.21 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.227 on 29581 degrees of freedom
## (229 observations deleted due to missingness)
## Multiple R-squared: 0.1656, Adjusted R-squared: 0.1655
## F-statistic: 2935 on 2 and 29581 DF, p-value: < 2.2e-16x1=seq(min(biov1$velocity, na.rm = T),max(biov1$velocity, na.rm = T),le=100)
p1=predict(lm0,newdata=data.frame(velocity=x1),type="response")
points(p1~x1,type="l",col=2,lwd=2) #bimodal relationship
10 All shifts vs same direction
ggplot(biov1, aes(x = velocity, y = SHIFT))+
geom_point(aes(color = lag2), alpha = .1)+
scale_color_viridis_c()+
geom_smooth(method = "lm")+
theme_classic()+
facet_wrap(Param~Type+ECO, scales = "free")## `geom_smooth()` using formula = 'y ~ x'
ggplot(biov1 %>%
filter(shift_vel_sign == "pospos" | shift_vel_sign == "negneg"),
aes(x = abs(velocity), y = abs(SHIFT)))+
geom_point(aes(color = lag2), alpha = .1)+
scale_color_viridis_c()+
geom_smooth(method = "lm")+
theme_classic()+
facet_wrap(Param~Type+ECO, scales = "free")## `geom_smooth()` using formula = 'y ~ x'
ggplot(biov1, aes(x = velocity, y = SHIFT_cor))+
geom_point(aes(color = lag2), alpha = .1)+
scale_color_viridis_c()+
geom_smooth(method = "lm")+
theme_classic()+
facet_wrap(Param~Type+ECO, scales = "free")## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 229 rows containing non-finite values (`stat_smooth()`).## Warning: Removed 229 rows containing missing values (`geom_point()`).
ggplot(biov1 %>%
filter(shift_vel_sign == "pospos" | shift_vel_sign == "negneg"),
aes(x = abs(velocity), y = abs(SHIFT_cor)))+
geom_point(aes(color = lag2), alpha = .1)+
scale_color_viridis_c()+
geom_smooth(method = "lm")+
theme_classic()+
facet_wrap(Param~Type+ECO, scales = "free")## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 126 rows containing non-finite values (`stat_smooth()`).## Warning: Removed 126 rows containing missing values (`geom_point()`).
11 Load genetic diversity data
- Fonseca et al. 2023 Evo Letters: Mitochondrial DNA (for vertebrates, arachnids, and insects. Chloroplast DNA and nucleotide diversity (for plants). (π) was calculated to describe patterns of intraspecific genetic diversity.
# Fonseca et al. 2023
gen_div <- fread(here::here("Data/Fonseca_etal_2023_EvoLetters.txt"))
gen_div$spp <- gen_div$Species
gen_div$spp <- gsub("-"," ",gen_div$spp)
gen_div$spp <- gsub("_"," ",gen_div$spp)
gen_div$spp <- Clean_Names(gen_div$spp,
return_gen_sps = TRUE)
##############################
# Remove species identified to the genus level or cf. species
all_sps <- unique(gen_div$spp)
length(all_sps)
spgen <- sapply(all_sps, function(x){
tmp. <- strsplit(x, " ")[[1]]
any(tmp. == "sp" | tmp. == "spp" | grepl("sp[.]",tmp.) | grepl("spp[.]",tmp.) | tmp. == "cf[.]" | tmp. == "cf" | length(tmp.) == 1 | length(tmp.) > 2)
})
spgen <- which(spgen)
spgen <- all_sps[-spgen]
all(spgen %in% all_sps)
all_sps <- all_sps[which(all_sps %in% spgen)]
length(all_sps) # N species after filtering species identified to the species level
# filter dataset
gen_div <- gen_div[which(gen_div$spp %in% all_sps),]
length(unique(gen_div$spp)) # 3521911.1 Fix names at the genetic data
# Species to find
mycols <- c("reported_name_fixed","scientificName","kingdom","phylum","class","order","family","db_code")
sps_accepted_names <- data.frame(matrix(ncol = length(mycols), nrow = length(unique(all_sps))))
names(sps_accepted_names) <- mycols
sps_accepted_names$reported_name_fixed <- unique(all_sps)
tofind_ <- sps_accepted_names[which(is.na(sps_accepted_names$scientificName)),]
tofind_ <- unique(tofind_$reported_name_fixed)
tofind <- data.frame(matrix(nrow = length(tofind_), ncol = 8))
names(tofind) = c("scientificName", "kingdom", "phylum", "class", "order", "family", "db", "db_code")
tofind <- data.frame(species = tofind_, tofind)
tofind <- tofind %>%
mutate(across(everything(), as.character))
# retrieve sp names
sp_names_found <- Find_Sci_Names(sp_name = tofind$species)
# ----------------
# Summary
# ----------------
#
# N taxa:
# 36697
# N taxa found:
# |db | N|
# |:----|-----:|
# |GBIF | 35903|
# |ITIS | 76|
# |NCBI | 278|
# N taxa not found:
# 440
## Add found species names to the sps_accepted_names
all(sp_names_found$requested_name %in% sps_accepted_names$reported_name_fixed)
for(i in 1:length(sp_names_found$requested_name)){
tofill <- unique(which(sps_accepted_names$reported_name_fixed == sp_names_found$requested_name[i]))
sps_accepted_names$scientificName[tofill] <- sp_names_found$scientificName[i]
sps_accepted_names$kingdom[tofill] <- sp_names_found$kingdom[i]
sps_accepted_names$phylum[tofill] <-sp_names_found$phylum[i]
sps_accepted_names$class[tofill] <- sp_names_found$class[i]
sps_accepted_names$order[tofill] <- sp_names_found$order[i]
sps_accepted_names$family[tofill] <- sp_names_found$family[i]
sps_accepted_names$db[tofill] <- sp_names_found$db[i]
sps_accepted_names$db_code[tofill] <- sp_names_found$db_code[i]
sps_accepted_names$method[tofill] <- sp_names_found$method[i]
}
sps_accepted_names$spp = sps_accepted_names$reported_name_fixed
## Keep original taxa information for those species we could not find a name at GBIF
pos <- which(is.na(sps_accepted_names$scientificName))
sps_accepted_names$scientificName[pos] <- sps_accepted_names$spp[pos]
for(i in 1:nrow(gen_div)){
gen_div$spp_new[i] <- sps_accepted_names$scientificName[which(sps_accepted_names$spp == gen_div$spp[i])]
}
write.csv(gen_div,
here::here('Data/Fonseca_etal_2023_EvoLetters_harmo.txt'), row.names = FALSE)12 Load harmonized genetic data
gen_div <- fread(here::here('Data/Fonseca_etal_2023_EvoLetters_harmo.txt'))12.1 Stats
# Gen
gen_sps <- unique(gen_div$spp_new)
bsi_v1 <- intersect(biov1$spp, gen_sps)
# Break down
tot1 <- data.frame(N_species_with_gen_data = length(gen_sps),
v1_matches = length(bsi_v1))
tot1## N_species_with_gen_data v1_matches
## 1 34702 2960# N_species_with_gen_data v1_matches
# 34702 2982
# Select only species in Bioshifts
gen_div <- gen_div %>% dplyr::filter(spp %in% bsi_v1)12.2 How many species with >1 gen div measurements ?
N_obs <- gen_div %>%
group_by(spp) %>%
tally() %>%
dplyr::filter(n>1)
# total N
length(unique(gen_div$spp))## [1] 2914# total N with > 1 observation
nrow(N_obs)## [1] 2# these species
DT::datatable(N_obs)# How many species with >1 gen div measurements at the same site?
# Check how many obs per site and marker type exist
N_obs <- gen_div %>%
group_by(spp, round(Longitude,1), round(Latitude,1)) %>%
tally() %>%
dplyr::filter(n>1)
DT::datatable(N_obs)No observations in the same site
12.3 Which species?
toplot_ <- gen_div %>%
group_by(Group) %>%
dplyr::summarise(Obs = length(spp),
Species = length(unique(spp))) %>%
gather(var, N, -c(Group))
ggplot(toplot_, aes(x = Group, y = N))+
geom_bar(stat="identity")+
geom_text(aes(y=N, label=N,
hjust = ifelse(N<max(N),-.1,1.1)), vjust=0.2, size=3,
position = position_dodge(0.9))+
theme_classic()+
coord_flip()+
facet_wrap(.~var, scales = "free", ncol=1)
12.4 Distribution of He
unique(gen_div$Group)## [1] "Actinopterygii" "Insecta" "Arachnida" "Magnoliopsida"
## [5] "Aves" "Mammalia" "Amphibia" "Reptilia"
## [9] "Pinopsida" "Polypodiopsida" "Liliopsida" "Bryopsida"
## [13] "Lycopodiopsida" "Elasmobranchii" "Psilotopsida" "Holocephali"
## [17] "Florideophyceae"to_plot <- gen_div %>%
select(Group,Nucleotide_diversity, TajimasD) %>%
filter(Group %in% c("Actinopterygii","Insecta","Aves", "Magnoliopsida", "Mammalia")) %>%
gather(metric, value, -Group)
ggplot(to_plot, aes(value, fill = Group, color = Group))+
geom_density(alpha = .3)+
scale_color_viridis_d()+
scale_fill_viridis_d()+
facet_wrap(.~metric, scales = "free")
12.5 High He with latitude?
12.5.1 Absolution latitude
toplot <- gen_div %>%
mutate(Hemisphere = ifelse(Latitude<0, "South", "North"),
lat_abs = abs(Latitude),
Nucleotide_diversity = scale(Nucleotide_diversity),
TajimasD = scale(TajimasD))
ggplot(toplot, aes(x = lat_abs, y = Nucleotide_diversity, color = Group))+
geom_point(alpha = .1)+
stat_smooth(method = "lm")+
facet_wrap(.~Group,scales = "free")## `geom_smooth()` using formula = 'y ~ x'## Warning in qt((1 - level)/2, df): NaNs produced## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf
mod1 <- lm(Nucleotide_diversity~lat_abs*Group, toplot)
emtrends(mod1, ~ lat_abs*Group, var = "lat_abs")## lat_abs Group lat_abs.trend SE df lower.CL upper.CL
## 47 Actinopterygii -0.02001 0.00540 2884 -0.0306 -0.00942
## 47 Amphibia -0.01105 0.01918 2884 -0.0487 0.02656
## 47 Arachnida -0.00856 0.01281 2884 -0.0337 0.01655
## 47 Aves -0.02014 0.00329 2884 -0.0266 -0.01369
## 47 Bryopsida 0.00515 0.03441 2884 -0.0623 0.07262
## 47 Elasmobranchii 0.00251 0.01710 2884 -0.0310 0.03605
## 47 Florideophyceae nonEst NA NA NA NA
## 47 Holocephali 0.00843 0.10358 2884 -0.1947 0.21151
## 47 Insecta -0.00825 0.00251 2884 -0.0132 -0.00333
## 47 Liliopsida 0.00274 0.00670 2884 -0.0104 0.01587
## 47 Lycopodiopsida -0.00808 0.08012 2884 -0.1652 0.14902
## 47 Magnoliopsida -0.00205 0.00522 2884 -0.0123 0.00819
## 47 Mammalia -0.02342 0.01687 2884 -0.0565 0.00965
## 47 Pinopsida 0.00341 0.04155 2884 -0.0781 0.08488
## 47 Polypodiopsida 0.00185 0.02034 2884 -0.0380 0.04174
## 47 Psilotopsida nonEst NA NA NA NA
## 47 Reptilia -0.03945 0.05755 2884 -0.1523 0.07339
##
## Confidence level used: 0.9512.6 Spatial pattern
mundi <- terra::vect(rnaturalearth::ne_coastline(scale = 110, returnclass = "sp"))## Warning: The `returnclass` argument of `ne_download()` sp as of rnaturalearth 1.0.0.
## ℹ Please use `sf` objects with {rnaturalearth}, support for Spatial objects
## (sp) will be removed in a future release of the package.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.# get vect gen div
vect_div <- terra::vect(gen_div %>%
dplyr::filter(spp %in% biov1$sp_name_std_v1),
geom=c("Longitude","Latitude"),crs=crs(mundi))
# get vect centroids bioshifts
vect_biov1 <- terra::vect(biov1 %>%
filter(spp %in% gen_div$spp),
geom=c("centroid.x","centroid.y"),crs=crs(mundi))
# get raster bioshifts shp files study areas
get_raster_bioshifts = "NO"
if(get_raster_bioshifts=="YES"){
my_ext = terra::ext(mundi)
my_crs = crs(mundi)
rast_biov1 <- terra::rast(my_ext, crs = my_crs, res = 0.5)
values(rast_biov1) <- 0
fgdb <- "C:/Users/brunn/ShadowDrive/CreateGeodatabaseBioShifts/Data/Study_Areas_v1/Study_Areas.gdb"
fc_list <- terra::vector_layers(fgdb)
fc_list <- fc_list[which(fc_list %in% unique(as.data.frame(vect_biov1)$ID))]
for(i in 1:length(fc_list)){ cat("\r",i,"from",length(fc_list))
tmp = terra::vect(fgdb, layer = fc_list[i])
tmp = terra::cells(rast_biov1,tmp)
tmp_cell = tmp[,2]
tmp_vals = rast_biov1[tmp_cell][,1]
rast_biov1[tmp_cell] = tmp_vals+1
}
names(rast_biov1) <- "SA"
rast_biov1[rast_biov1==0] <- NA
writeRaster(rast_biov1, here::here("Data/raster_bioshifts_SA.tif"), overwrite = TRUE)
} else {
rast_biov1 <- terra::rast(here::here("Data/raster_bioshifts_SA.tif"))
}
values_biov1 <- na.omit(values(rast_biov1))
ggplot()+
ggtitle("Genetic diversity data")+
geom_spatvector(data=mundi)+
geom_spatvector(data=vect_div)+
theme_blank()
ggplot()+
ggtitle("Genetic diversity data + Bioshifts study areas")+
geom_spatraster(data=rast_biov1)+
scale_fill_whitebox_b(
palette = "muted",
na.value = "white",
breaks = seq(min(values_biov1), max(values_biov1), 1))+
geom_spatvector(data=vect_div, alpha = .2)+
geom_spatvector(data=mundi)+
theme_blank()
13 Merge genetic and bioshifts data
#merge all data frames in list
gen_data_v1_avg <- append(list(gen_div[,c(3:9,15)]), list(biov1)) %>% reduce(full_join, by='spp')
gen_data_v1_avg <- gen_data_v1_avg %>% dplyr::filter(!is.na(Nucleotide_diversity), !is.na(SHIFT))
# N species
length(unique(gen_data_v1_avg$spp))## [1] 291414 Check if gen data overlap with study areas
fgdb <- "C:/Users/brunn/ShadowDrive/CreateGeodatabaseBioShifts/Data/Study_Areas_v1/Study_Areas.gdb"
SAs <- terra::vector_layers(fgdb)
SAs <- SAs[which(SAs %in% unique(as.data.frame(vect_biov1)$ID))]
overlap_SA <- data.frame()
for(i in 1:length(SAs)){ cat("\r", i ,"from", length(SAs))
SA_i <- SAs[i]
# get species in v1 within SA i
v1_sps_in_SA_i <- unique(gen_data_v1_avg$spp[which(gen_data_v1_avg$ID == SA_i)])
# get the coordinates of genetic data for those species
gen_coord <- gen_data_v1_avg %>%
dplyr::filter(spp %in% v1_sps_in_SA_i) %>%
dplyr::select(c("spp","Longitude","Latitude")) %>%
dplyr::filter(!duplicated(spp)) %>%
vect(geom = c("Longitude","Latitude"))
# load SA i
SA_i_shp <- terra::vect(fgdb, layer = SA_i)
# which one of these species overlap with the SA?
tmp = c(relate(ext(SA_i_shp), gen_coord, "contains"))
SA_i <- data.frame(spp = gen_coord$spp,
ID = SA_i,
overlap_SA = ifelse(tmp, 1, 0))
overlap_SA <- rbind(overlap_SA,
SA_i)
}##
1 from 203## Warning: [vect] Z coordinates ignored##
2 from 203## Warning: [vect] Z coordinates ignored##
3 from 203## Warning: [vect] Z coordinates ignored##
4 from 203## Warning: [vect] Z coordinates ignored##
5 from 203## Warning: [vect] Z coordinates ignored##
6 from 203## Warning: [vect] Z coordinates ignored##
7 from 203## Warning: [vect] Z coordinates ignored##
8 from 203
9 from 203## Warning: [vect] Z coordinates ignored##
10 from 203## Warning: [vect] Z coordinates ignored##
11 from 203## Warning: [vect] Z coordinates ignored##
12 from 203
13 from 203
14 from 203
15 from 203
16 from 203
17 from 203
18 from 203
19 from 203
20 from 203
21 from 203## Warning: [vect] Z coordinates ignored##
22 from 203## Warning: [vect] Z coordinates ignored##
23 from 203## Warning: [vect] Z coordinates ignored##
24 from 203
25 from 203
26 from 203
27 from 203
28 from 203
29 from 203
30 from 203
31 from 203## Warning: [vect] Z coordinates ignored##
32 from 203
33 from 203
34 from 203## Warning: [vect] Z coordinates ignored##
35 from 203## Warning: [vect] Z coordinates ignored##
36 from 203## Warning: [vect] Z coordinates ignored##
37 from 203## Warning: [vect] Z coordinates ignored##
38 from 203## Warning: [vect] Z coordinates ignored##
39 from 203## Warning: [vect] Z coordinates ignored##
40 from 203## Warning: [vect] Z coordinates ignored##
41 from 203## Warning: [vect] Z coordinates ignored##
42 from 203## Warning: [vect] Z coordinates ignored##
43 from 203## Warning: [vect] Z coordinates ignored##
44 from 203## Warning: [vect] Z coordinates ignored##
45 from 203## Warning: [vect] Z coordinates ignored##
46 from 203## Warning: [vect] Z coordinates ignored##
47 from 203## Warning: [vect] Z coordinates ignored##
48 from 203## Warning: [vect] Z coordinates ignored##
49 from 203## Warning: [vect] Z coordinates ignored##
50 from 203## Warning: [vect] Z coordinates ignored##
51 from 203## Warning: [vect] Z coordinates ignored##
52 from 203
53 from 203
54 from 203
55 from 203
56 from 203
57 from 203
58 from 203
59 from 203## Warning: [vect] Z coordinates ignored##
60 from 203## Warning: [vect] Z coordinates ignored##
61 from 203## Warning: [vect] Z coordinates ignored##
62 from 203## Warning: [vect] Z coordinates ignored##
63 from 203
64 from 203## Warning: [vect] Z coordinates ignored##
65 from 203
66 from 203
67 from 203
68 from 203
69 from 203## Warning: [vect] Z coordinates ignored##
70 from 203## Warning: [vect] Z coordinates ignored##
71 from 203## Warning: [vect] Z coordinates ignored##
72 from 203## Warning: [vect] Z coordinates ignored##
73 from 203## Warning: [vect] Z coordinates ignored##
74 from 203
75 from 203
76 from 203## Warning: [vect] Z coordinates ignored##
77 from 203
78 from 203
79 from 203
80 from 203
81 from 203## Warning: [vect] Z coordinates ignored##
82 from 203## Warning: [vect] Z coordinates ignored##
83 from 203## Warning: [vect] Z coordinates ignored##
84 from 203## Warning: [vect] Z coordinates ignored##
85 from 203
86 from 203## Warning: [vect] Z coordinates ignored##
87 from 203## Warning: [vect] Z coordinates ignored##
88 from 203## Warning: [vect] Z coordinates ignored##
89 from 203## Warning: [vect] Z coordinates ignored##
90 from 203## Warning: [vect] Z coordinates ignored##
91 from 203
92 from 203## Warning: [vect] Z coordinates ignored##
93 from 203## Warning: [vect] Z coordinates ignored##
94 from 203
95 from 203
96 from 203
97 from 203
98 from 203## Warning: [vect] Z coordinates ignored##
99 from 203## Warning: [vect] Z coordinates ignored##
100 from 203## Warning: [vect] Z coordinates ignored##
101 from 203## Warning: [vect] Z coordinates ignored##
102 from 203## Warning: [vect] Z coordinates ignored##
103 from 203## Warning: [vect] Z coordinates ignored##
104 from 203## Warning: [vect] Z coordinates ignored##
105 from 203
106 from 203
107 from 203
108 from 203
109 from 203## Warning: [vect] Z coordinates ignored##
110 from 203
111 from 203
112 from 203## Warning: [vect] Z coordinates ignored##
113 from 203## Warning: [vect] Z coordinates ignored##
114 from 203## Warning: [vect] Z coordinates ignored##
115 from 203## Warning: [vect] Z coordinates ignored##
116 from 203
117 from 203## Warning: [vect] Z coordinates ignored##
118 from 203## Warning: [vect] Z coordinates ignored##
119 from 203## Warning: [vect] Z coordinates ignored##
120 from 203
121 from 203
122 from 203## Warning: [vect] Z coordinates ignored##
123 from 203
124 from 203## Warning: [vect] Z coordinates ignored##
125 from 203## Warning: [vect] Z coordinates ignored##
126 from 203## Warning: [vect] Z coordinates ignored##
127 from 203
128 from 203## Warning: [vect] Z coordinates ignored##
129 from 203## Warning: [vect] Z coordinates ignored##
130 from 203## Warning: [vect] Z coordinates ignored##
131 from 203## Warning: [vect] Z coordinates ignored##
132 from 203## Warning: [vect] Z coordinates ignored##
133 from 203## Warning: [vect] Z coordinates ignored##
134 from 203## Warning: [vect] Z coordinates ignored##
135 from 203## Warning: [vect] Z coordinates ignored##
136 from 203## Warning: [vect] Z coordinates ignored##
137 from 203## Warning: [vect] Z coordinates ignored##
138 from 203## Warning: [vect] Z coordinates ignored##
139 from 203## Warning: [vect] Z coordinates ignored##
140 from 203## Warning: [vect] Z coordinates ignored##
141 from 203## Warning: [vect] Z coordinates ignored##
142 from 203## Warning: [vect] Z coordinates ignored##
143 from 203## Warning: [vect] Z coordinates ignored##
144 from 203## Warning: [vect] Z coordinates ignored##
145 from 203## Warning: [vect] Z coordinates ignored##
146 from 203## Warning: [vect] Z coordinates ignored##
147 from 203## Warning: [vect] Z coordinates ignored##
148 from 203## Warning: [vect] Z coordinates ignored##
149 from 203## Warning: [vect] Z coordinates ignored##
150 from 203## Warning: [vect] Z coordinates ignored##
151 from 203## Warning: [vect] Z coordinates ignored##
152 from 203## Warning: [vect] Z coordinates ignored##
153 from 203## Warning: [vect] Z coordinates ignored##
154 from 203## Warning: [vect] Z coordinates ignored##
155 from 203## Warning: [vect] Z coordinates ignored##
156 from 203## Warning: [vect] Z coordinates ignored##
157 from 203## Warning: [vect] Z coordinates ignored##
158 from 203## Warning: [vect] Z coordinates ignored##
159 from 203## Warning: [vect] Z coordinates ignored##
160 from 203
161 from 203## Warning: [vect] Z coordinates ignored##
162 from 203## Warning: [vect] Z coordinates ignored##
163 from 203## Warning: [vect] Z coordinates ignored##
164 from 203## Warning: [vect] Z coordinates ignored##
165 from 203## Warning: [vect] Z coordinates ignored##
166 from 203## Warning: [vect] Z coordinates ignored##
167 from 203## Warning: [vect] Z coordinates ignored##
168 from 203## Warning: [vect] Z coordinates ignored##
169 from 203## Warning: [vect] Z coordinates ignored##
170 from 203## Warning: [vect] Z coordinates ignored##
171 from 203## Warning: [vect] Z coordinates ignored##
172 from 203## Warning: [vect] Z coordinates ignored##
173 from 203## Warning: [vect] Z coordinates ignored##
174 from 203## Warning: [vect] Z coordinates ignored##
175 from 203## Warning: [vect] Z coordinates ignored##
176 from 203## Warning: [vect] Z coordinates ignored##
177 from 203## Warning: [vect] Z coordinates ignored##
178 from 203## Warning: [vect] Z coordinates ignored##
179 from 203
180 from 203## Warning: [vect] Z coordinates ignored##
181 from 203## Warning: [vect] Z coordinates ignored##
182 from 203## Warning: [vect] Z coordinates ignored##
183 from 203## Warning: [vect] Z coordinates ignored##
184 from 203## Warning: [vect] Z coordinates ignored##
185 from 203## Warning: [vect] Z coordinates ignored##
186 from 203## Warning: [vect] Z coordinates ignored##
187 from 203## Warning: [vect] Z coordinates ignored##
188 from 203## Warning: [vect] Z coordinates ignored##
189 from 203## Warning: [vect] Z coordinates ignored##
190 from 203## Warning: [vect] Z coordinates ignored##
191 from 203## Warning: [vect] Z coordinates ignored##
192 from 203## Warning: [vect] Z coordinates ignored##
193 from 203## Warning: [vect] Z coordinates ignored##
194 from 203## Warning: [vect] Z coordinates ignored##
195 from 203## Warning: [vect] Z coordinates ignored##
196 from 203## Warning: [vect] Z coordinates ignored##
197 from 203## Warning: [vect] Z coordinates ignored##
198 from 203## Warning: [vect] Z coordinates ignored##
199 from 203## Warning: [vect] Z coordinates ignored##
200 from 203## Warning: [vect] Z coordinates ignored##
201 from 203## Warning: [vect] Z coordinates ignored##
202 from 203## Warning: [vect] Z coordinates ignored##
203 from 203## Warning: [vect] Z coordinates ignoredtable(overlap_SA$overlap_SA)##
## 0 1
## 7698 714gen_data_v1_avg <- merge(gen_data_v1_avg,
overlap_SA,
by=c("spp","ID"))14.1 Save data
write.csv(gen_data_v1_avg,here::here("Data/gen_data_final_m2_fonseca.csv"),row.names = FALSE)15 Explore data
# n range shifts
nrow(gen_data_v1_avg)## [1] 1070915.1 Which species?
# N species
length(unique(gen_data_v1_avg$spp))## [1] 2914# Use only plus plus or minus minus
table(gen_data_v1_avg$shift_vel_sign)##
## negneg negpos posneg pospos
## 191 3915 284 6319# N species
length(unique(gen_data_v1_avg$spp))## [1] 2914toplot_ <- gen_data_v1_avg %>%
group_by(Group) %>%
dplyr::summarise(Species = length(unique(spp)),
Obs = length(spp))
toplot_ <- reshape2::melt(toplot_)## Using Group as id variablesggplot(toplot_, aes(x = Group, y = value))+
geom_bar(stat="identity")+
geom_text(aes(y=value, label=value),
hjust = -.1, vjust=0.2, size=3,
position = position_dodge(0.9))+
theme_classic()+
coord_flip()+
facet_wrap(.~variable, scales = "free", ncol = 1)
15.2 Compare lag values of the full bioshifts vs the subset (genetic) dataset
toplot <- rbind(
data.frame(dataset = "Full",
lags = biov1$lag2,
shifts = biov1$SHIFT),
data.frame(dataset = "Subset",
lags = gen_data_v1_avg$lag2,
shifts = gen_data_v1_avg$SHIFT))
# lags
ggplot(toplot, aes(x = dataset, y=lags, fill = dataset, color = dataset))+
geom_boxplot(alpha = .5, outlier.shape = NA)+
scale_y_continuous(limits = quantile(toplot$lags, c(0.1, 0.9), na.rm = T))## Warning: Removed 7893 rows containing non-finite values (`stat_boxplot()`).
tapply(toplot$lags, toplot$dataset, summary)## $Full
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -145.1653 -0.4375 0.6237 0.3779 2.1868 13.1259
##
## $Subset
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -145.16531 -0.97750 0.51824 0.09254 2.15833 13.12588mod1 <- aov(lags~dataset, data = toplot)
summary(mod1)## Df Sum Sq Mean Sq F value Pr(>F)
## dataset 1 642 641.7 28.16 1.12e-07 ***
## Residuals 40520 923357 22.8
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# shift
ggplot(toplot, aes(x = dataset, y=shifts, fill = dataset, color = dataset))+
geom_boxplot(alpha = .5, outlier.shape = NA)+
scale_y_continuous(limits = quantile(toplot$shifts, c(0.1, 0.9), na.rm = T))## Warning: Removed 8100 rows containing non-finite values (`stat_boxplot()`).
tapply(toplot$shifts, toplot$dataset, summary)## $Full
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -119.1696 -0.4468 0.3056 1.1671 2.4345 146.3000
##
## $Subset
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -94.6628 -0.3000 0.3867 1.3493 2.7400 146.3000mod1 <- aov(shifts~dataset, data = toplot)
summary(mod1)## Df Sum Sq Mean Sq F value Pr(>F)
## dataset 1 262 261.58 8.673 0.00323 **
## Residuals 40520 1222023 30.16
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1There are difference in mean values of lags, but no difference in mean shift values, between the full and subset datasets.
15.3 Lag
15.3.1 General pattern
gen_data_v1_avg %>%
dplyr::filter(Type == "LAT") %>%
ggplot(aes(x = Nucleotide_diversity, y = lagC))+
ggtitle("Latitude")+
xlab("Genetic diversity")+
ylab("lagC (Velocity - Shift)") +
scale_color_viridis_c(direction = -1)+
geom_point(alpha = .5, size = 3, aes(color = SHIFT_abs))+
theme_bw()+
geom_smooth(method = "lm", color = "black", size=1)+
stat_smooth(geom="line", method = "lm", alpha=0.3, se = FALSE,
aes(by=as.factor(Article_ID)), color = "black")## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning in stat_smooth(geom = "line", method = "lm", alpha = 0.3, se = FALSE, :
## Ignoring unknown aesthetics: by## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 46 rows containing non-finite values (`stat_smooth()`).## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 46 rows containing non-finite values (`stat_smooth()`).## Warning: Removed 46 rows containing missing values (`geom_point()`).
gen_data_v1_avg %>%
dplyr::filter(Type == "ELE") %>%
ggplot(aes(x = Nucleotide_diversity, y = lagC, color = SHIFT_abs))+
ggtitle("Elevation")+
xlab("Genetic diversity")+
ylab("lagC (Velocity - Shift)") +
scale_color_viridis_c(direction = -1)+
geom_point(alpha = .5, size = 3)+
theme_bw()+
geom_smooth(method = "lm", color = "black", size=1)+
stat_smooth(geom="line", method = "lm", alpha=0.3, se = FALSE,
aes(by=as.factor(Article_ID)), color = "black")## Warning in stat_smooth(geom = "line", method = "lm", alpha = 0.3, se = FALSE, :
## Ignoring unknown aesthetics: by## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 1 rows containing non-finite values (`stat_smooth()`).## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 1 rows containing non-finite values (`stat_smooth()`).## Warning: Removed 1 rows containing missing values (`geom_point()`).
15.3.2 By Param
15.3.2.1 Nucleotide_diversity
gen_data_v1_avg %>%
dplyr::filter(Type == "LAT") %>%
ggplot(aes(x = Nucleotide_diversity, y = lagC, color = SHIFT_abs))+
ggtitle("Latitude")+
xlab("Genetic diversity")+
ylab("lagC (Velocity - Shift)") +
scale_color_viridis_c(direction = -1)+
geom_point(alpha = .5, size = 3)+
theme_bw()+
geom_smooth(method = "lm", color = "black", size=1)+
stat_smooth(geom="line", method = "lm", alpha=0.3, se = FALSE,
aes(by=as.factor(Article_ID)), color = "black")+
facet_wrap(.~Param, scales = "free", ncol = 3)## Warning in stat_smooth(geom = "line", method = "lm", alpha = 0.3, se = FALSE, :
## Ignoring unknown aesthetics: by## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 46 rows containing non-finite values (`stat_smooth()`).## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 46 rows containing non-finite values (`stat_smooth()`).## Warning: Removed 46 rows containing missing values (`geom_point()`).
gen_data_v1_avg %>%
dplyr::filter(Type == "ELE") %>%
ggplot(aes(x = Nucleotide_diversity, y = lagC, color = SHIFT_abs))+
ggtitle("Elevation")+
xlab("Genetic diversity")+
ylab("lagC (Velocity - Shift)") +
scale_color_viridis_c(direction = -1)+
geom_point(alpha = .5, size = 3)+
theme_bw()+
geom_smooth(method = "lm", color = "black", size=1)+
stat_smooth(geom="line", method = "lm", alpha=0.3, se = FALSE,
aes(by=as.factor(Article_ID)), color = "black")+
facet_wrap(.~Param, scales = "free", ncol = 3)## Warning in stat_smooth(geom = "line", method = "lm", alpha = 0.3, se = FALSE, :
## Ignoring unknown aesthetics: by## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 1 rows containing non-finite values (`stat_smooth()`).## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 1 rows containing non-finite values (`stat_smooth()`).## Warning: Removed 1 rows containing missing values (`geom_point()`).
15.3.3 By Group
gen_data_v1_avg %>%
dplyr::filter(Type == "LAT") %>%
ggplot(aes(x = Nucleotide_diversity, y = lagC, color = SHIFT_abs))+
ggtitle("Latitude")+
xlab("Genetic diversity")+
ylab("lagC (Velocity - Shift)") +
scale_color_viridis_c(direction = -1)+
geom_point(alpha = .5, size = 3)+
theme_bw()+
geom_smooth(method = "lm", color = "black", size=1)+
stat_smooth(geom="line", method = "lm", alpha=0.3, se = FALSE,
aes(by=as.factor(Article_ID)), color = "black")+
facet_wrap(.~Group, scales = "free", ncol = 3)## Warning in stat_smooth(geom = "line", method = "lm", alpha = 0.3, se = FALSE, :
## Ignoring unknown aesthetics: by## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 46 rows containing non-finite values (`stat_smooth()`).## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 46 rows containing non-finite values (`stat_smooth()`).## Warning: Removed 46 rows containing missing values (`geom_point()`).
gen_data_v1_avg %>%
dplyr::filter(Type == "ELE") %>%
ggplot(aes(x = Nucleotide_diversity, y = lagC, color = SHIFT_abs))+
ggtitle("Elevation")+
xlab("Genetic diversity")+
ylab("lagC (Velocity - Shift)") +
scale_color_viridis_c(direction = -1)+
geom_point(alpha = .5, size = 3)+
theme_bw()+
geom_smooth(method = "lm", color = "black", size=1)+
stat_smooth(geom="line", method = "lm", alpha=0.3, se = FALSE,
aes(by=as.factor(Article_ID)), color = "black")+
facet_wrap(.~Group, scales = "free", ncol = 3)## Warning in stat_smooth(geom = "line", method = "lm", alpha = 0.3, se = FALSE, :
## Ignoring unknown aesthetics: by## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 1 rows containing non-finite values (`stat_smooth()`).## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 1 rows containing non-finite values (`stat_smooth()`).## Warning: Removed 1 rows containing missing values (`geom_point()`).
15.3.4 By Group / Param
gen_data_v1_avg %>%
dplyr::filter(Type == "LAT") %>%
ggplot(aes(x = Nucleotide_diversity, y = lagC, color = SHIFT_abs))+
ggtitle("Latitude")+
xlab("Genetic diversity")+
ylab("lagC (Velocity - Shift)") +
scale_color_viridis_c(direction = -1)+
geom_point(alpha = .5, size = 3)+
theme_bw()+
geom_smooth(method = "lm", color = "black", size=1)+
stat_smooth(geom="line", method = "lm", alpha=0.3, se = FALSE,
aes(by=as.factor(Article_ID)), color = "black")+
facet_grid(Group~Param,scales = "free")## Warning in stat_smooth(geom = "line", method = "lm", alpha = 0.3, se = FALSE, :
## Ignoring unknown aesthetics: by## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 46 rows containing non-finite values (`stat_smooth()`).## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 46 rows containing non-finite values (`stat_smooth()`).## Warning: Removed 46 rows containing missing values (`geom_point()`).
gen_data_v1_avg %>%
dplyr::filter(Type == "ELE") %>%
ggplot(aes(x = Nucleotide_diversity, y = lagC, color = SHIFT_abs))+
ggtitle("Elevation")+
xlab("Genetic diversity")+
ylab("lagC (Velocity - Shift)") +
scale_color_viridis_c(direction = -1)+
geom_point(alpha = .5, size = 3)+
theme_bw()+
geom_smooth(method = "lm", color = "black", size=1)+
stat_smooth(geom="line", method = "lm", alpha=0.3, se = FALSE,
aes(by=as.factor(Article_ID)), color = "black")+
facet_grid(Group~Param,scales = "free")## Warning in stat_smooth(geom = "line", method = "lm", alpha = 0.3, se = FALSE, :
## Ignoring unknown aesthetics: by## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 1 rows containing non-finite values (`stat_smooth()`).## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 1 rows containing non-finite values (`stat_smooth()`).## Warning: Removed 1 rows containing missing values (`geom_point()`).
15.4 Lag 2
15.4.1 General pattern
gen_data_v1_avg %>%
dplyr::filter(Type == "LAT") %>%
ggplot(aes(x = Nucleotide_diversity, y = lagC2, color = SHIFT_abs))+
ggtitle("Latitude")+
xlab("Genetic diversity")+
ylab("lagC (Velocity - Shift)") +
scale_color_viridis_c(direction = -1)+
geom_point(alpha = .5, size = 3)+
theme_bw()+
geom_smooth(method = "lm", color = "black", size=1)+
stat_smooth(geom="line", method = "lm", alpha=0.3, se = FALSE,
aes(by=as.factor(Article_ID)), color = "black")## Warning in stat_smooth(geom = "line", method = "lm", alpha = 0.3, se = FALSE, :
## Ignoring unknown aesthetics: by## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 46 rows containing non-finite values (`stat_smooth()`).## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 46 rows containing non-finite values (`stat_smooth()`).## Warning: Removed 46 rows containing missing values (`geom_point()`).
gen_data_v1_avg %>%
dplyr::filter(Type == "ELE") %>%
ggplot(aes(x = Nucleotide_diversity, y = lagC2, color = SHIFT_abs))+
ggtitle("Elevation")+
xlab("Genetic diversity")+
ylab("lagC (Velocity - Shift)") +
scale_color_viridis_c(direction = -1)+
geom_point(alpha = .5, size = 3)+
theme_bw()+
geom_smooth(method = "lm", color = "black", size=1)+
stat_smooth(geom="line", method = "lm", alpha=0.3, se = FALSE,
aes(by=as.factor(Article_ID)), color = "black")## Warning in stat_smooth(geom = "line", method = "lm", alpha = 0.3, se = FALSE, :
## Ignoring unknown aesthetics: by## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 1 rows containing non-finite values (`stat_smooth()`).## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 1 rows containing non-finite values (`stat_smooth()`).## Warning: Removed 1 rows containing missing values (`geom_point()`).
15.4.2 By Param
gen_data_v1_avg %>%
dplyr::filter(Type == "LAT") %>%
ggplot(aes(x = Nucleotide_diversity, y = lagC2, color = SHIFT_abs))+
ggtitle("Latitude")+
xlab("Genetic diversity")+
ylab("lagC (Velocity - Shift)") +
scale_color_viridis_c(direction = -1)+
geom_point(alpha = .5, size = 3)+
theme_bw()+
geom_smooth(method = "lm", color = "black", size=1)+
stat_smooth(geom="line", method = "lm", alpha=0.3, se = FALSE,
aes(by=as.factor(Article_ID)), color = "black")+
facet_wrap(.~Param, scales = "free", ncol = 3)## Warning in stat_smooth(geom = "line", method = "lm", alpha = 0.3, se = FALSE, :
## Ignoring unknown aesthetics: by## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 46 rows containing non-finite values (`stat_smooth()`).## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 46 rows containing non-finite values (`stat_smooth()`).## Warning: Removed 46 rows containing missing values (`geom_point()`).
gen_data_v1_avg %>%
dplyr::filter(Type == "ELE") %>%
ggplot(aes(x = Nucleotide_diversity, y = lagC2, color = SHIFT_abs))+
ggtitle("Elevation")+
xlab("Genetic diversity")+
ylab("lagC (Velocity - Shift)") +
scale_color_viridis_c(direction = -1)+
geom_point(alpha = .5, size = 3)+
theme_bw()+
geom_smooth(method = "lm", color = "black", size=1)+
stat_smooth(geom="line", method = "lm", alpha=0.3, se = FALSE,
aes(by=as.factor(Article_ID)), color = "black")+
facet_wrap(.~Param, scales = "free", ncol = 3)## Warning in stat_smooth(geom = "line", method = "lm", alpha = 0.3, se = FALSE, :
## Ignoring unknown aesthetics: by## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 1 rows containing non-finite values (`stat_smooth()`).## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 1 rows containing non-finite values (`stat_smooth()`).## Warning: Removed 1 rows containing missing values (`geom_point()`).
15.4.3 By Group
gen_data_v1_avg %>%
dplyr::filter(Type == "LAT") %>%
ggplot(aes(x = Nucleotide_diversity, y = lagC2, color = SHIFT_abs))+
ggtitle("Latitude")+
xlab("Genetic diversity")+
ylab("lagC (Velocity - Shift)") +
scale_color_viridis_c(direction = -1)+
geom_point(alpha = .5, size = 3)+
theme_bw()+
geom_smooth(method = "lm", color = "black", size=1)+
stat_smooth(geom="line", method = "lm", alpha=0.3, se = FALSE,
aes(by=as.factor(Article_ID)), color = "black")+
facet_wrap(.~Group, scales = "free", ncol = 3)## Warning in stat_smooth(geom = "line", method = "lm", alpha = 0.3, se = FALSE, :
## Ignoring unknown aesthetics: by## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 46 rows containing non-finite values (`stat_smooth()`).## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 46 rows containing non-finite values (`stat_smooth()`).## Warning: Removed 46 rows containing missing values (`geom_point()`).
gen_data_v1_avg %>%
dplyr::filter(Type == "ELE") %>%
ggplot(aes(x = Nucleotide_diversity, y = lagC2, color = SHIFT_abs))+
ggtitle("Elevation")+
xlab("Genetic diversity")+
ylab("lagC (Velocity - Shift)") +
scale_color_viridis_c(direction = -1)+
geom_point(alpha = .5, size = 3)+
theme_bw()+
geom_smooth(method = "lm", color = "black", size=1)+
stat_smooth(geom="line", method = "lm", alpha=0.3, se = FALSE,
aes(by=as.factor(Article_ID)), color = "black")+
facet_wrap(.~Group, scales = "free", ncol = 3)## Warning in stat_smooth(geom = "line", method = "lm", alpha = 0.3, se = FALSE, :
## Ignoring unknown aesthetics: by## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 1 rows containing non-finite values (`stat_smooth()`).## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 1 rows containing non-finite values (`stat_smooth()`).## Warning: Removed 1 rows containing missing values (`geom_point()`).
15.4.4 By Group / Param
gen_data_v1_avg %>%
dplyr::filter(Type == "LAT") %>%
ggplot(aes(x = Nucleotide_diversity, y = lagC2, color = SHIFT_abs))+
ggtitle("Latitude")+
xlab("Genetic diversity")+
ylab("lagC (Velocity - Shift)") +
scale_color_viridis_c(direction = -1)+
geom_point(alpha = .5, size = 3)+
theme_bw()+
geom_smooth(method = "lm", color = "black", size=1)+
stat_smooth(geom="line", method = "lm", alpha=0.3, se = FALSE,
aes(by=as.factor(Article_ID)), color = "black")+
facet_wrap(Param~Group, scales = "free")## Warning in stat_smooth(geom = "line", method = "lm", alpha = 0.3, se = FALSE, :
## Ignoring unknown aesthetics: by## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 46 rows containing non-finite values (`stat_smooth()`).## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 46 rows containing non-finite values (`stat_smooth()`).## Warning: Removed 46 rows containing missing values (`geom_point()`).
gen_data_v1_avg %>%
dplyr::filter(Type == "ELE") %>%
ggplot(aes(x = Nucleotide_diversity, y = lagC2, color = SHIFT_abs))+
ggtitle("Elevation")+
xlab("Genetic diversity")+
ylab("lagC (Velocity - Shift)") +
scale_color_viridis_c(direction = -1)+
geom_point(alpha = .5, size = 3)+
theme_bw()+
geom_smooth(method = "lm", color = "black", size=1)+
stat_smooth(geom="line", method = "lm", alpha=0.3, se = FALSE,
aes(by=as.factor(Article_ID)), color = "black")+
facet_wrap(Param~Group, scales = "free")## Warning in stat_smooth(geom = "line", method = "lm", alpha = 0.3, se = FALSE, :
## Ignoring unknown aesthetics: by## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 1 rows containing non-finite values (`stat_smooth()`).## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 1 rows containing non-finite values (`stat_smooth()`).## Warning: Removed 1 rows containing missing values (`geom_point()`).
15.5 SHIFT
15.5.1 General pattern
# gen_data_v1_avg <- gen_data_v1_avg_fon
# gen_data_v1_avg$Nucleotide_diversity <- gen_data_v1_avg$He
gen_data_v1_avg %>%
dplyr::filter(Type == "LAT") %>%
ggplot(aes(x = Nucleotide_diversity, y = sqrt(abs(SHIFT_cor)), color = lag))+
ggtitle("Latitude")+
xlab("Genetic diversity")+
ylab("Range shift (km/year)") +
scale_color_viridis_c(direction = -1)+
geom_point(alpha = .5, size = 3)+
theme_bw()+
geom_smooth(method = "lm", color = "black", size=1)+
stat_smooth(geom="line", method = "lm", alpha=0.3, se = FALSE,
aes(by=as.factor(Article_ID)), color = "black")## Warning in stat_smooth(geom = "line", method = "lm", alpha = 0.3, se = FALSE, :
## Ignoring unknown aesthetics: by## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 46 rows containing non-finite values (`stat_smooth()`).## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 46 rows containing non-finite values (`stat_smooth()`).## Warning: Removed 46 rows containing missing values (`geom_point()`).
gen_data_v1_avg %>%
dplyr::filter(Type == "ELE") %>%
ggplot(aes(x = Nucleotide_diversity, y = sqrt(abs(SHIFT_cor)), color = lag))+
ggtitle("Elevation")+
xlab("Genetic diversity")+
ylab("Range shift (km/year)") +
scale_color_viridis_c(direction = -1)+
geom_point(alpha = .5, size = 3)+
theme_bw()+
geom_smooth(method = "lm", color = "black", size=1)+
stat_smooth(geom="line", method = "lm", alpha=0.3, se = FALSE,
aes(by=as.factor(Article_ID)), color = "black")## Warning in stat_smooth(geom = "line", method = "lm", alpha = 0.3, se = FALSE, :
## Ignoring unknown aesthetics: by## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 1 rows containing non-finite values (`stat_smooth()`).## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 1 rows containing non-finite values (`stat_smooth()`).## Warning: Removed 1 rows containing missing values (`geom_point()`).
15.5.2 By Param
gen_data_v1_avg %>%
dplyr::filter(Type == "LAT") %>%
ggplot(aes(x = Nucleotide_diversity, y = sqrt(abs(SHIFT_cor)), color = lag))+
ggtitle("Latitude")+
xlab("Genetic diversity")+
ylab("Range shift (km/year)") +
scale_color_viridis_c(direction = -1)+
geom_point(alpha = .5, size = 3)+
theme_bw()+
geom_smooth(method = "lm", color = "black", size=1)+
stat_smooth(geom="line", method = "lm", alpha=0.3, se = FALSE,
aes(by=as.factor(Article_ID)), color = "black")+
facet_wrap(.~Param, scales = "free", ncol = 3)## Warning in stat_smooth(geom = "line", method = "lm", alpha = 0.3, se = FALSE, :
## Ignoring unknown aesthetics: by## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 46 rows containing non-finite values (`stat_smooth()`).## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 46 rows containing non-finite values (`stat_smooth()`).## Warning: Removed 46 rows containing missing values (`geom_point()`).
gen_data_v1_avg %>%
dplyr::filter(Type == "ELE") %>%
ggplot(aes(x = Nucleotide_diversity, y = sqrt(abs(SHIFT_cor)), color = lag))+
ggtitle("Elevation")+
xlab("Genetic diversity")+
ylab("Range shift (km/year)") +
scale_color_viridis_c(direction = -1)+
geom_point(alpha = .5, size = 3)+
theme_bw()+
geom_smooth(method = "lm", color = "black", size=1)+
stat_smooth(geom="line", method = "lm", alpha=0.3, se = FALSE,
aes(by=as.factor(Article_ID)), color = "black")+
facet_wrap(.~Param, scales = "free", ncol = 3)## Warning in stat_smooth(geom = "line", method = "lm", alpha = 0.3, se = FALSE, :
## Ignoring unknown aesthetics: by## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 1 rows containing non-finite values (`stat_smooth()`).## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 1 rows containing non-finite values (`stat_smooth()`).## Warning: Removed 1 rows containing missing values (`geom_point()`).
15.5.3 By Group
gen_data_v1_avg %>%
dplyr::filter(Type == "LAT") %>%
ggplot(aes(x = Nucleotide_diversity, y = sqrt(abs(SHIFT_cor)), color = lag))+
ggtitle("Latitude")+
xlab("Genetic diversity")+
ylab("Range shift (km/year)") +
scale_color_viridis_c(direction = -1)+
geom_point(alpha = .5, size = 3)+
theme_bw()+
geom_smooth(method = "lm", color = "black", size=1)+
stat_smooth(geom="line", method = "lm", alpha=0.3, se = FALSE,
aes(by=as.factor(Article_ID)), color = "black")+
facet_wrap(.~Group, scales = "free", ncol = 3)## Warning in stat_smooth(geom = "line", method = "lm", alpha = 0.3, se = FALSE, :
## Ignoring unknown aesthetics: by## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 46 rows containing non-finite values (`stat_smooth()`).## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 46 rows containing non-finite values (`stat_smooth()`).## Warning: Removed 46 rows containing missing values (`geom_point()`).
gen_data_v1_avg %>%
dplyr::filter(Type == "ELE") %>%
ggplot(aes(x = Nucleotide_diversity, y = sqrt(abs(SHIFT_cor)), color = lag))+
ggtitle("Elevation")+
xlab("Genetic diversity")+
ylab("Range shift (km/year)") +
scale_color_viridis_c(direction = -1)+
geom_point(alpha = .5, size = 3)+
theme_bw()+
geom_smooth(method = "lm", color = "black", size=1)+
stat_smooth(geom="line", method = "lm", alpha=0.3, se = FALSE,
aes(by=as.factor(Article_ID)), color = "black")+
facet_wrap(.~Group, scales = "free", ncol = 3)## Warning in stat_smooth(geom = "line", method = "lm", alpha = 0.3, se = FALSE, :
## Ignoring unknown aesthetics: by## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 1 rows containing non-finite values (`stat_smooth()`).## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 1 rows containing non-finite values (`stat_smooth()`).## Warning: Removed 1 rows containing missing values (`geom_point()`).
15.5.4 By Group / Param
gen_data_v1_avg %>%
dplyr::filter(Type == "LAT") %>%
ggplot(aes(x = Nucleotide_diversity, y = (abs(SHIFT_cor)), color = lag))+
ggtitle("Latitude")+
xlab("Genetic diversity")+
ylab("Range shift (km/year)") +
scale_color_viridis_c(direction = -1)+
geom_point(alpha = .5, size = 3)+
theme_bw()+
geom_smooth(method = "lm", color = "black", size=1)+
stat_smooth(geom="line", method = "lm", alpha=0.3, se = FALSE,
aes(by=as.factor(Article_ID)), color = "black")+
facet_wrap(Param~Group, scales = "free")## Warning in stat_smooth(geom = "line", method = "lm", alpha = 0.3, se = FALSE, :
## Ignoring unknown aesthetics: by## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 46 rows containing non-finite values (`stat_smooth()`).## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 46 rows containing non-finite values (`stat_smooth()`).## Warning: Removed 46 rows containing missing values (`geom_point()`).
gen_data_v1_avg %>%
dplyr::filter(Type == "ELE") %>%
ggplot(aes(x = Nucleotide_diversity, y = (abs(SHIFT_cor)), color = lag))+
ggtitle("Elevation")+
xlab("Genetic diversity")+
ylab("Range shift (km/year)") +
scale_color_viridis_c(direction = -1)+
geom_point(alpha = .5, size = 3)+
theme_bw()+
geom_smooth(method = "lm", color = "black", size=1)+
stat_smooth(geom="line", method = "lm", alpha=0.3, se = FALSE,
aes(by=as.factor(Article_ID)), color = "black")+
facet_wrap(Param~Group, scales = "free")## Warning in stat_smooth(geom = "line", method = "lm", alpha = 0.3, se = FALSE, :
## Ignoring unknown aesthetics: by## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 1 rows containing non-finite values (`stat_smooth()`).## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 1 rows containing non-finite values (`stat_smooth()`).## Warning: Removed 1 rows containing missing values (`geom_point()`).